#include "redmath.h"

extern "C"
void
load_module ()
{
  addToSystem ("pi", new dNumber (3.14159265));
  addToSystem ("gcd", dGCD);
  addToSystem ("exp", dExp);
  addToSystem ("pow", dPow);
}

#ifndef BIGNUM
/* The Binary GCD algorithm, taken from en.wikipedia.org/wiki/Binary_GCD_algorithm */
int gcd(int u, int v)
{
    int shift;

    /* GCD(0,x) := x */
    if (u == 0 || v == 0)
      return u | v;

    /* Let shift := lg K, where K is the greatest power of 2
       dividing both u and v. */
    for (shift = 0; ((u | v) & 1) == 0; ++shift) {
        u >>= 1;
        v >>= 1;
    }

    while ((u & 1) == 0)
      u >>= 1;

    /* From here on, u is always odd. */
    do {
        while ((v & 1) == 0)
          v >>= 1;

        /* Now u and v are both odd, so diff(u, v) is even.
           Let u = min(u, v), v = diff(u, v)/2. */
        if (u < v) {
            v -= u;
        } else {
            int diff = u - v;
            u = v;
            v = diff;
        }
        v >>= 1;
    } while (v != 0);

    return u << shift;
}

#endif

dObject *
dGCD (dList *args)
{
  if (args == 0 || args->size != 2
      || (*args)[0]->type != 'i'
      || (*args)[1]->type != 'i')
    return new dError ("gcd requires two integers");
#ifdef BIGNUM
  mpz_class result, x, y;
  x = ((dNumber*)(*args)[0])->Integer ();
  y = ((dNumber*)(*args)[1])->Integer ();
  mpz_gcd (result.get_mpz_t (), x.get_mpz_t (), y.get_mpz_t ());
  return new dNumber (result);
#else
  int x, y;
  x = ((dNumber*)(*args)[0])->Integer ();
  y = ((dNumber*)(*args)[1])->Integer ();
  return new dNumber (gcd (x, y));
#endif
}

dObject *dExp (dList *args)
{
  if (args == 0 || args->size != 1)
    return new dError ("exp requires one number");
#ifdef BIGNUM
  mpf_class x = ((dNumber*)(*args)[0])->Decimal ();
  return new dNumber (exp (x.get_d ()));
#else
  double x = ((dNumber*)(*args)[0])->Decimal ();
  return new dNumber (exp (x));
#endif
}

dObject *dPow (dList *args)
{
  if (args == 0 || args->size != 2)
    return new dError ("pow requires two numbers");
#ifdef BIGNUM
  mpf_class result, x;
  unsigned long int y;
  x = ((dNumber*)(*args)[0])->Decimal ();
  y = ((dNumber*)(*args)[1])->Decimal ().get_si ();
  mpf_pow_ui (result.get_mpf_t (), x.get_mpf_t (), y);
  return new dNumber (result);
#else
  double x, y;
  x = ((dNumber*)(*args)[0])->Decimal ();
  y = ((dNumber*)(*args)[1])->Decimal ();
  return new dNumber (pow (x, y));
#endif
}
